%     VIRTUAL HC - программа для получения эффективных свойств сотовых заполнителей
%     Copyright (C) 2018-2020  Михаил Козлов
%     Сайт проекта: <http://www.fastenerdoctor.ru/>
%     Адрес электронной почты для связи с разработчиком:
%     <doctor.fastener@gmail.com>
% 
%     Эта программа является свободным программным обеспечением: Вы можете
%     распространять ее и/или изменять, соблюдая условия Универсальной общественной
%     лицензии GNU, опубликованной Фондом свободного программного обеспечения,
%     либо редакции 3, либо (на Ваше усмотрение) любой редакции, выпущенной позже.
% 
%     Эта программа распространяется в надежде на то, что она окажется полезной,
%     но БЕЗ КАКИХ-ЛИБО ГАРАНТИЙ, включая даже подразумеваемую гарантию ТОВАРНОЙ ПРИГОДНОСТИ,
%     либо ПРИГОДНОСТИ ДЛЯ ОПРЕДЕЛЕННЫХ ЦЕЛЕЙ. Ознакомьтесь с Универсальной общественной
%     лицензией GNU для получения более подробной информации.
% 
%     Вы должны были получить копию Универсальной общественной лицензии GNU
%     вместе с этой программой. Если Вы ее не получили, то перейдите по адресу:
%     <http://www.gnu.org/licenses/>.
% 
%     This program is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
% 
%     This program is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
%     GNU General Public License for more details.
% 
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <https://www.gnu.org/licenses/>.

function [NODES, ELEMENTS, RVE_Coordinates, Angle, Node_Merge_Tolerance] = Generate_Mesh(RVE_Type, RVE_Frame_Points, RVE_Size, Cell_Mesh_Type, Element_Order, N_Elements_Edge, Element_Critical_Ratio, Foil_Thickness, Cell_Filler, RVE_Array_L, RVE_Array_W, RVE_Array_T, Analysis_Type, RVE_High_Model, N_Elements_High)

global ErrCode;

global Merge_Tolerance_Multiplier;

NODES = [];
ELEMENTS = [];
RVE_Coordinates = [];
Node_Merge_Tolerance = [];

fprintf (1, 'Построение сетки...                                     ');

% Вычисление угла раскрытия сторон ячейки
switch Cell_Mesh_Type
    
    case 'SHELL'
        
        switch RVE_Type
            
            case 'TYPE_A'
                
                Vec_1 = RVE_Frame_Points(7,:) - RVE_Frame_Points(6,:);
                Vec_2 = RVE_Frame_Points(10,:) - RVE_Frame_Points(6,:);
                
                Angle = acosd(dot(Vec_1, Vec_2) / (norm(Vec_1) * norm(Vec_2)));
                
            case 'TYPE_B'
                
                Vec_1 = RVE_Frame_Points(11,:) - RVE_Frame_Points(10,:);
                Vec_2 = RVE_Frame_Points(2,:) - RVE_Frame_Points(10,:);
                
                Angle = acosd(dot(Vec_1, Vec_2) / (norm(Vec_1) * norm(Vec_2)));
                
            otherwise
                
                fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое значение параметра RVE_Type !\n');
                ErrCode = -1;
                return;
            
        end
        
    case 'SOLID'
        
        switch RVE_Type
            
            case 'TYPE_A'
                
                Vec_1 = RVE_Frame_Points(15,:) - RVE_Frame_Points(13,:);
                Vec_2 = RVE_Frame_Points(22,:) - RVE_Frame_Points(13,:);
                
                Angle = acosd(dot(Vec_1, Vec_2) / (norm(Vec_1) * norm(Vec_2)));
                
            case 'TYPE_B'
                
                Vec_1 = RVE_Frame_Points(3,:) - RVE_Frame_Points(2,:);
                Vec_2 = RVE_Frame_Points(10,:) - RVE_Frame_Points(2,:);
                
                Angle = acosd(dot(Vec_1, Vec_2) / (norm(Vec_1) * norm(Vec_2)));
                
            otherwise
                
                fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое значение параметра RVE_Type !\n');
                ErrCode = -1;
                return;
            
        end
        
    otherwise
        
        fprintf (1, '\nОшибка функции Generate_Mesh: недопустимое значение параметра Cell_Mesh_Type !\n');
        ErrCode = -1;
        return;
    
end

% Выбор алгоритма для построения плоской сетки
if (Angle < 45)
    
    [NODES_2D, ELEMENTS_2D, Min_Element_Edge, N_Layers_Suggestion] = Mesher_2D_Type_1(RVE_Type, RVE_Frame_Points, Cell_Mesh_Type, N_Elements_Edge, Element_Critical_Ratio, Foil_Thickness, Cell_Filler);

    if (ErrCode < 0)
        return;
    end
    
end

if ((Angle >= 45) && (Angle <= 90))
    
    [NODES_2D, ELEMENTS_2D, Min_Element_Edge, N_Layers_Suggestion] = Mesher_2D_Type_2(RVE_Type, RVE_Frame_Points, Cell_Mesh_Type, N_Elements_Edge, Element_Critical_Ratio, Foil_Thickness, Cell_Filler);

    if (ErrCode < 0)
        return;
    end
    
end

if ((Angle > 90) && (Angle <= 135))
    
    [NODES_2D, ELEMENTS_2D, Min_Element_Edge, N_Layers_Suggestion] = Mesher_2D_Type_3(RVE_Type, RVE_Frame_Points, Cell_Mesh_Type, N_Elements_Edge, Element_Critical_Ratio, Foil_Thickness, Cell_Filler);

    if (ErrCode < 0)
        return;
    end
    
end

if (Angle > 135)
    
    [NODES_2D, ELEMENTS_2D, Min_Element_Edge, N_Layers_Suggestion] = Mesher_2D_Type_4(RVE_Type, RVE_Frame_Points, Cell_Mesh_Type, N_Elements_Edge, Element_Critical_Ratio, Foil_Thickness, Cell_Filler);

    if (ErrCode < 0)
        return;
    end
    
end

% Число элементов по высоте ячейки в случае автоматически вычисленной высоты
if (strcmp(RVE_High_Model, 'AUTO'))
    N_Elements_High = N_Layers_Suggestion;
end

% Увеличение числа элементов по высоте ячейки в случае превышения максимально допустимого соотношения сторон элемента
Element_Ratio = RVE_Size(3)/N_Layers_Suggestion/Min_Element_Edge;

if (Element_Ratio > Element_Critical_Ratio)
    N_Elements_High = ceil(N_Elements_High * Element_Ratio / Element_Critical_Ratio);
end

% В механическом расчете должно быть не менее двух элементов по высоте, иначе не будет узлов для некоторых сэтов
if (strcmp(Analysis_Type, 'MECHANICAL') && (N_Elements_High == 1))
    N_Elements_High = 2;
end

% % Четное число элементов по высоте в тепловом расчете необходимо для создания сэта RVE_CenterTemp (он не нужен для осреднения)
% if (strcmp(Analysis_Type, 'HEAT_TRANSFER') && (rem(N_Elements_High,2) == 1))
%     N_Elements_High = N_Elements_High + 1;
% end

% Обновление минимальной длины ребра элемента
Min_Element_Edge = min(Min_Element_Edge, RVE_Size(3) / N_Elements_High);

% Допуск для слиянии узлов
Node_Merge_Tolerance = Min_Element_Edge * Merge_Tolerance_Multiplier;

% Число узлов плоской сетки
N_Nodes_2D = size(NODES_2D,1);

% Выделение памяти под узлы объемной сетки
NODES = zeros((N_Elements_High + 1) * N_Nodes_2D, 4);

% Заполнение массива узлов объемной сетки
for i=0:N_Elements_High
    
    NODES(i*N_Nodes_2D + 1:(i+1)*N_Nodes_2D,:) = NODES_2D(:,:);
    
    if (i > 0)
        
        NODES(i*N_Nodes_2D + 1:(i+1)*N_Nodes_2D, 1) = NODES((i-1)*N_Nodes_2D + 1:i*N_Nodes_2D, 1) + ones(N_Nodes_2D,1) * N_Nodes_2D;
        NODES(i*N_Nodes_2D + 1:(i+1)*N_Nodes_2D, 4) = ones(N_Nodes_2D,1) * i/N_Elements_High * RVE_Size(3);
        
    end
    
end

% Число элементов плоской сетки
N_Elements_2D = size(ELEMENTS_2D,1);

% Выделение памяти элементы объемной сетки
ELEMENTS = zeros(N_Elements_High * N_Elements_2D, 12);

% Заполнение первого слоя объемной сетки
for i=1:N_Elements_2D
    
    switch ELEMENTS_2D(i,2)
        
        case 1
            
            ELEMENTS(i,1) = ELEMENTS_2D(i,1);
            ELEMENTS(i,2) = 2;
            ELEMENTS(i,3) = 4;
            ELEMENTS(i,4) = ELEMENTS_2D(i,4);
            
            ELEMENTS(i,5) = ELEMENTS_2D(i,6);
            ELEMENTS(i,6) = ELEMENTS_2D(i,5);
            ELEMENTS(i,7) = ELEMENTS_2D(i,5) + N_Nodes_2D;
            ELEMENTS(i,8) = ELEMENTS_2D(i,6) + N_Nodes_2D;
            
        case 2
            
            switch ELEMENTS_2D(i,3)
                
                case 3
                    
                    ELEMENTS(i,1) = ELEMENTS_2D(i,1);
                    ELEMENTS(i,2) = 3;
                    ELEMENTS(i,3) = 6;
                    ELEMENTS(i,4) = ELEMENTS_2D(i,4);
                    
                    Element_Normal = cross(NODES_2D(ELEMENTS_2D(i,6),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4), NODES_2D(ELEMENTS_2D(i,7),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4));
                    
                    if (dot(Element_Normal, [0,0,1]) > 0)
                        
                        ELEMENTS(i,5) = ELEMENTS_2D(i,5);
                        ELEMENTS(i,6) = ELEMENTS_2D(i,6);
                        ELEMENTS(i,7) = ELEMENTS_2D(i,7);
                        ELEMENTS(i,8) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                        ELEMENTS(i,9) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                        ELEMENTS(i,10) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                        
                    else
                        
                        ELEMENTS(i,5) = ELEMENTS_2D(i,7);
                        ELEMENTS(i,6) = ELEMENTS_2D(i,6);
                        ELEMENTS(i,7) = ELEMENTS_2D(i,5);
                        ELEMENTS(i,8) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                        ELEMENTS(i,9) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                        ELEMENTS(i,10) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                        
                    end
                    
                case 4
                    
                    ELEMENTS(i,1) = ELEMENTS_2D(i,1);
                    ELEMENTS(i,2) = 3;
                    ELEMENTS(i,3) = 8;
                    ELEMENTS(i,4) = ELEMENTS_2D(i,4);
                    
                    Element_Normal = cross(NODES_2D(ELEMENTS_2D(i,6),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4), NODES_2D(ELEMENTS_2D(i,8),2:4) - NODES_2D(ELEMENTS_2D(i,5),2:4));

                    if (dot(Element_Normal, [0,0,1]) > 0)
                        
                        ELEMENTS(i,5) = ELEMENTS_2D(i,5);
                        ELEMENTS(i,6) = ELEMENTS_2D(i,6);
                        ELEMENTS(i,7) = ELEMENTS_2D(i,7);
                        ELEMENTS(i,8) = ELEMENTS_2D(i,8);
                        ELEMENTS(i,9) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                        ELEMENTS(i,10) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                        ELEMENTS(i,11) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                        ELEMENTS(i,12) = ELEMENTS_2D(i,8) + N_Nodes_2D;
                        
                    else
                        
                        ELEMENTS(i,5) = ELEMENTS_2D(i,8);
                        ELEMENTS(i,6) = ELEMENTS_2D(i,7);
                        ELEMENTS(i,7) = ELEMENTS_2D(i,6);
                        ELEMENTS(i,8) = ELEMENTS_2D(i,5);
                        ELEMENTS(i,9) = ELEMENTS_2D(i,8) + N_Nodes_2D;
                        ELEMENTS(i,10) = ELEMENTS_2D(i,7) + N_Nodes_2D;
                        ELEMENTS(i,11) = ELEMENTS_2D(i,6) + N_Nodes_2D;
                        ELEMENTS(i,12) = ELEMENTS_2D(i,5) + N_Nodes_2D;
                        
                    end
                
            end
        
    end
    
end

% Заполнение остальных слоев объемной сетки
Mask = ELEMENTS(1 : N_Elements_2D, 5:12) > 0;
Mask = double(Mask);
    
for i=1:N_Elements_High-1
    
    ELEMENTS(i*N_Elements_2D + 1 : (i+1)*N_Elements_2D, 1) = ELEMENTS((i-1)*N_Elements_2D + 1 : i*N_Elements_2D, 1) + ones(N_Elements_2D,1)*N_Elements_2D;
    ELEMENTS(i*N_Elements_2D + 1 : (i+1)*N_Elements_2D, 2:4) = ELEMENTS((i-1)*N_Elements_2D + 1 : i*N_Elements_2D, 2:4);
    ELEMENTS(i*N_Elements_2D + 1 : (i+1)*N_Elements_2D, 5:12) = ELEMENTS((i-1)*N_Elements_2D + 1 : i*N_Elements_2D, 5:12) + Mask*N_Nodes_2D;
    
end

% Габариты ячейки
RVE_Coordinates = zeros(3,2);

RVE_Coordinates(1,1) = -RVE_Size(1)/2;
RVE_Coordinates(1,2) = RVE_Size(1)/2;
RVE_Coordinates(2,1) = -RVE_Size(2)/2;
RVE_Coordinates(2,2) = RVE_Size(2)/2;
RVE_Coordinates(3,1) = 0;
RVE_Coordinates(3,2) = RVE_Size(3);

fprintf (1, 'OK!\n');

% Конвертация элементов в элементы второго порядка
if (strcmp(Element_Order, 'QUADRATIC'))
    
    fprintf (1, 'Конвертация элементов в квадратичные...                 ');
    
    % Уменьшение допуска для слияния узлов
    Node_Merge_Tolerance = Node_Merge_Tolerance/2;
    
    [NODES, ELEMENTS] = Convert_Mesh_To_Quad(NODES, ELEMENTS, Node_Merge_Tolerance);
    
    fprintf (1, 'OK!\n');
    
end

% Создание массива ячеек
if ((RVE_Array_L > 1) || (RVE_Array_W > 1) || (RVE_Array_T > 1))
    
    fprintf (1, 'Создание массива ячеек...                               ');
    
    [NODES, ELEMENTS, RVE_Coordinates] = Create_RVE_Array(NODES, ELEMENTS, Node_Merge_Tolerance, RVE_Array_L, RVE_Array_W, RVE_Array_T, RVE_Size, RVE_Coordinates, Cell_Mesh_Type);
    
    fprintf (1, 'OK!\n');
    
end
